      !<@file 
      !>@brief PDMDFDF1: Problema Direto, Método das Diferenças Finitas, Domínio da Frequência, 1: Hustedt et al, 2004, equação 8.
      !>
      !> Baseado em Hustedt, B.; Operto, S. & Virieux, J. (2004), 'Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling', Geophysical Journal International 157(3), 1269-1296.
      subroutine PDMDFDF1(omega,iflagF,iflagG,iflagH)


      use globais
      implicit none

      double precision :: omega

      complex*16 :: c1, c2, c3, c4, c5 ! variaveis auxiliares

      integer :: iflagF
      integer :: iflagG
      integer :: iflagH
      
      integer :: kk,noL,noK

      integer::knnz,iL,jL


!       para o solver do lapack ou superlu
      integer :: NRHS, info

! !	para solver lapack
!       complex, dimension(:,:), allocatable :: LwB ! matriz do sistema em banda
!       integer :: KU,KL,LDAB
!       integer, dimension(:), allocatable :: IPIV !auxiliar para o solver
!       complex, dimension(:,:), allocatable :: f !lado direito do sistema 




! 	para a SuperLU
!       complex*16 ::f(maxNeq,maxNsource)
!       integer :: rowind(maxNnz), colptr(maxNeq+1) !CCS      
!       complex*16 ::  values(maxNnz)
      integer ::  ldf, iopt
      integer*8 :: factors




      integer :: linha, col1,col2,col3,col4,col5 !  variáveis auxiliares para a montagem da matriz do sistema
      integer :: ndonoh, col !  variáveis auxiliares para prescrever superficie livre
      integer :: i,j,cont

      complex*16 :: aux




      ! vetor csix e csiz, jah que estes variam para cada frequencia omega

      do i=1,Nx

	  if(i.le.Npml )then ! dentro da primeira camada de absorcao vertical
	      gammax = Cpml*cos(Pi/2.0 * dfloat(i-1)/dfloat(Npml-1))
	  elseif(i.gt.(Nx-Npml))then ! dentro da segunda camada de absorcao vertical
	      gammax = Cpml*cos(Pi/2.0 * dfloat(Nx-i)/dfloat(Npml-1))
	  else ! dentro do modelo
	      gammax = 0.0
	  endif

	  csix(i) = 1.0  + imag * gammax/omega

      enddo





      do j=1,Nz

! 	  if(j.le.Npml )then ! dentro da primeira camada de absorcao horizontal    (para duas camadas de absorcao horizontal)
! 	      gammaz = Cpml*cos(Pi/2.0 * dfloat(j-1)/dfloat(Npml-1))
! ! 	      gammaz = 0.0
! 	  elseif(j.gt.(Nz-Npml))then ! dentro da segunda camada de absorcao horizontal
! 	      gammaz = Cpml*cos(Pi/2.0 * dfloat(Nz-j)/dfloat(Npml-1))
! 	  else ! dentro do modelo
! 	      gammaz = 0.0
! 	  endif

	  if(j.gt.(Nz-Npml))then ! dentro da segunda camada de absorcao horizontal
	      gammaz = Cpml*cos(Pi/2.0 * dfloat(Nz-j)/dfloat(Npml-1))
	  else ! dentro do modelo
	      gammaz = 0.0
	  endif

	  csiz(j) = 1.0  + imag * gammaz/omega
	
      enddo

!       open(unit=2,file='perfil.txt')
!       do j=1,Nz
! 	  write(2,*)c(Nx/2,j)
!       enddo
!       close(2)

      Neq = Nx*Nz
      if(Neq.gt.maxNeq)then	
	  write(*,*)'Neq.gt.maxNeq',Neq,'>',maxNeq
      endif

! 	!para lapack
! 	KL = Nx
! 	KU = Nx
! 	LDAB = 2*KL+KU+1
! 
! 	write(*,*)'alocando LwB...'
! 	allocate(LwB(LDAB,Neq))
! 
! 	write(*,*)'alocando f...'
! 	allocate(f(Neq,Nsource))






!       "staggered-grid five-point stencil" 
! 	pags 1273: mixed-grid and staggered-grid finite-difference methods for 
!       frequency-domain acoustic wave modelling (hustedt, operto, virieux, 2004)






      write(*,*)'#Montando matrix do sistema...',kw,'/',Nomega,omega/2/pi,' hz'


      nnz=0 ! numero de nao zeros


      j=1  ! condicao de contorno , nohs da superficie
      do i=1, Nx

	  linha = (j-1)*Nx + i   !numero do noh ij
	  col1 = (j-1)*Nx + i	   ! corresponde ao noh i,j  !numero do noh ij
! 	    Lw(linha,col1) = 1.0d0
! 	    LwB(KL+KU+1+linha-col1,col1)=1.0d0    !coloca 1 no elemento da diagonal
  
	  !estou armazenando a matriz transposta com ccs
	  nnz= nnz+1 
	  colptr(linha)=nnz ! guardo o nnz do primeiro nao zero de cada linha
	  rowind(nnz)=col1 ! indice da linha eh o numero da col
	  values(nnz) = 1.0d0 !valor nao nulo
      enddo
	    

      do j=2, Nz

	  do i=1, Nx

	      linha = (j-1)*Nx + i   !numero do noh ij

	      ! ordem das colunas, da menor p maior: col4, col2, col1, col3, col5
	      col1 = (j-1)*Nx + i	   ! corresponde ao noh i,j  !numero do noh ij
	      col2 = (j-1)*Nx + (i-1)  ! corresponde ao noh i-1,j !numero do noh da esquerda
	      col3 = (j-1)*Nx + (i+1)  ! corresponde ao noh i+1,j	!numero do noh da direita
	      col4 = ((j-1)-1)*Nx + i ! corresponde ao noh i,j-1  !numero do noh de cima
	      col5 = ((j+1)-1)*Nx + i  ! corresponde ao noh i,j+1  !numero do noh de baixo
! 		write(*,*)'#   ',i,j,col1,col2,col3,col4,col5

	      if(i.ne.1 .and. i.ne.Nx .and.j.ne.1 .and. j.ne.Nz )then
		  ! staggered grid five point stencil (hustedt et al., 2004 equacao 8)
		  c2 = (1.0d0/(csix(i) * delta**2))  * (b(i,j)+b(i-1,j) )/(csix(i)+csix(i-1))    
		  c3 = (1.0d0/(csix(i) * delta**2))  * (b(i,j)+b(i+1,j) )/(csix(i)+csix(i+1))  
		  c4 = (1.0d0/(csiz(j) * delta**2))  * (b(i,j)+b(i,j-1) )/(csiz(j)+csiz(j-1))
		  c5 = (1.0d0/(csiz(j) * delta**2))  * (b(i,j)+b(i,j+1) )/(csiz(j)+csiz(j+1)) 
		  c1 = omega**2/K(i,j) - c2 -c3 - c4 - c5


! 		  ! Diferencas finitas de segunda ordem sem pml
! 		  c2=1.0d0/(delta**2 * rho(i-1,j))
! 		  c3=1.0d0/(delta**2 * rho(i+1,j))
! 		  c4=1.0d0/(delta**2 * rho(i,j-1))
! 		  c5=1.0d0/(delta**2 * rho(i,j+1)) 
! 		  c1 = omega**2/(rho(i,j)*c(i,j)**2) - c2 -c3 - c4 - c5

! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col2) = c2
! 		  Lw(linha,col3) = c3
! 		  Lw(linha,col4) = c4
! 		  Lw(linha,col5) = c5
! 		  LwB(KL+KU+1+linha-col1,col1) = c1
! 		  LwB(KL+KU+1+linha-col2,col2) = c2
! 		  LwB(KL+KU+1+linha-col3,col3) = c3
! 		  LwB(KL+KU+1+linha-col4,col4) = c4
! 		  LwB(KL+KU+1+linha-col5,col5) = c5

! 			!estou armazenando a matriz transposta com ccs
		  nnz= nnz+1 
		  colptr(linha)=nnz ! guardo o nnz do primeiro nao zero de cada linha
		  rowind(nnz)=col4 ! indice da linha eh o numero da col
		  values(nnz) = c4 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col2 ! indice da linha eh o numero da col
		  values(nnz) = c2 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col1 ! indice da linha eh o numero da col
		  values(nnz) = c1 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col3 ! indice da linha eh o numero da col
		  values(nnz) = c3 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col5 ! indice da linha eh o numero da col
		  values(nnz) = c5 !valor nao nulo

	      elseif(i.eq.1)then
		  c1 = -(1.0/delta + imag * omega/c(i,j))
		  c3 = +(1.0/delta)
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col3) = c3
! 		  LwB(KL+KU+1+linha-col1,col1) = c1
! 		  LwB(KL+KU+1+linha-col3,col3) = c3

! 			!estou armazenando a matriz transposta com ccs
		  nnz= nnz+1 
		  colptr(linha)=nnz ! guardo o nnz do primeiro nao zero de cada linha
		  rowind(nnz)=col1 ! indice da linha eh o numero da col
		  values(nnz) = c1 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col3 ! indice da linha eh o numero da col
		  values(nnz) = c3 !valor nao nulo


	      elseif(i.eq.Nx)then
		  c1 = +(1.0/delta + imag * omega/c(i,j))
		  c2 = -(1.0/delta)
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col2) = c2
! 		  LwB(KL+KU+1+linha-col1,col1) = c1
! 		  LwB(KL+KU+1+linha-col2,col2) = c2

! 			!estou armazenando a matriz transposta com ccs
		  nnz= nnz+1 
		  colptr(linha)=nnz ! guardo o nnz do primeiro nao zero de cada linha
		  rowind(nnz)=col2 ! indice da linha eh o numero da col
		  values(nnz) = c2 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col1 ! indice da linha eh o numero da col
		  values(nnz) = c1 !valor nao nulo


	      elseif(j.eq.1)then
		  c1 = -(1.0/delta + imag * omega/c(i,j))
		  c5 = +(1.0/delta)
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col5) = c5
! 		  LwB(KL+KU+1+linha-col1,col1) = c1
! 		  LwB(KL+KU+1+linha-col5,col5) = c5

! 			!estou armazenando a matriz transposta com ccs
		  nnz= nnz+1 
		  colptr(linha)=nnz ! guardo o nnz do primeiro nao zero de cada linha
		  rowind(nnz)=col1 ! indice da linha eh o numero da col
		  values(nnz) = c1 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col5 ! indice da linha eh o numero da col
		  values(nnz) = c5 !valor nao nulo

	      elseif(j.eq.Nz)then
		  c1 = +(1.0/delta + imag * omega/c(i,j))
		  c4 = -(1.0/delta)
! 		  Lw(linha,col1) = c1
! 		  Lw(linha,col4) = c4
! 		  LwB(KL+KU+1+linha-col1,col1) = c1
! 		  LwB(KL+KU+1+linha-col4,col4) = c4

! 			!estou armazenando a matriz transposta com ccs 
		  nnz= nnz+1 
		  colptr(linha)=nnz ! guardo o nnz do primeiro nao zero de cada linha
		  rowind(nnz)=col4 ! indice da linha eh o numero da col
		  values(nnz) = c4 !valor nao nulo
		  nnz= nnz+1 
		  rowind(nnz)=col1 ! indice da linha eh o numero da col
		  values(nnz) = c1 !valor nao nulo
	      endif
      
	  enddo
      enddo
      colptr(Neq+1)=nnz+1








! ! nao precisa mais prescrever deste jeito, depois da montagem da matriz... no loop de cima já se faz isso
! 	! para prescrever u=0 na superficie livre (j=1), substituir a linha e a coluna referente ao noh ij inteira por zero, exceto o elemento da diagonal que é feito = a 1, zerar no vetor independente a linha de cada no da superficie livre  
! 	j=1 !nohs da superficie livre
! 	do i=1,Nx
! 
! 	    j=1 !nohs da superficie livre
! 	    ndonoh  = (j-1)*Nx + i    ! numero do noh i,j
! 
! 	    linha = ndonoh
! 	    do col=1,Neq  ! zerando todos os elementos da linha "ndonoh"
! 		if((KL+KU+1+linha-col).ge.1 .and. (KL+KU+1+linha-col).le.LDAB .and. col .ge. 1 .and. col .le. Neq)then
! 		    LwB(KL+KU+1+linha-col,col)=0.0d0  
! 		endif
! 	    enddo
! 
! 	    col=ndonoh
! 	    do linha=1,Neq  ! zerando todos os elementos da coluna "ndonoh"
! 		if((KL+KU+1+linha-col).ge.1 .and. (KL+KU+1+linha-col).le.LDAB .and. col .ge. 1 .and. col .le. Neq)then
! 		    LwB(KL+KU+1+linha-col,col)=0.0d0  
! 		endif	
! 	    enddo
! 
! 	    linha = ndonoh	
! 	    col=ndonoh
! 	    LwB(KL+KU+1+linha-col,col)=1.0d0    !coloca 1 no elemento da diagonal	    
! 	    do ks=1,Nsource
! 		f(linha,ks)=0.0d0  !coloca 0 no vetor independente (nem precisa isso, ja q a fonte nunca estah em j=1, ou seja , ja vai ser zero)
! 	    enddo
! 	enddo
! 	!fim da condicao de contorno superficie livre





      if(iflagF.eq.1)then ! inicializar com 0 Misfit antes de entrar no loop dos grupos de pares fonte receptor
	  Misfit = 0.0d0
      endif
      if(iflagG.eq.1)then ! inicializar gradmisfit e matriz de precondicionamento
	  do noL=1,Nx0*Nz0 
	      gradmisfit(noL) = 0.0d0
	  enddo
      endif

      if(iflagH.eq.1)then
! 
! 	  if(ioptH.ne.0)then !< isto é, se já houve fatoração ou solução do sistema, deve-se limpar a memória alocada na fatoracao, depois atualizaer H
! 	      nrhsH = 1
! 	      ldyH = Nx0*nz0
! 	      ioptH = 3    !< Last, free the storage allocated inside SuperLU
! 	      call c_fortran_dgssv( ioptH, ldyH, Hnnz, nrhsH, Hvalue, Hrowind, Hcolptr, &
! 			  & f, ldyH, factorsH, infoH ) !CCS !aqui o f nao é usado, esta rotina nao conhece o vetor independente do sistema, coloquuei qq um
! 	      if (infoH .eq. 0) then
! 		  write (*,*) 'H free the storage succeeded'
! 	      else
! 		  write(*,*) 'H INFO from free the storage = ', infoH
! 		  stop
! 	      endif
! 	  endif



! 	  !matriz de precondicionamento CHEIA  !ccs, monta por coluna
! 	  Hnnz = 0 ! numero de elementos nao nulos da matriz de precondicionamento 
! 	  do jL=1,Nz0  !loop no grid de incognitas do pi
! 	      do iL=1,Nx0 
! 		  noL = (jL-1)*Nx0 + il !numero da incognita (coluna da Hessiana)
! 
! 		  noK = 1
! 		      Hnnz=Hnnz+1  !H(K,L)
! 		      Hcolptr(noL)= Hnnz	   ! nnz do primeiro nao nulo da coluna ! CCS da H 
! 		      Hrowind(Hnnz)= noK  ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL     ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0   ! CCS da H 
! 		  do noK = 2,Nx0*Nz0
! 		      Hnnz=Hnnz+1  !H(K,L)
! 		      Hrowind(Hnnz)= noK  ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL     ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0   ! CCS da H 
! 		  enddo
! 	      enddo
! 	  enddo
! 	  Hcolptr(Nx0*Nz0+1)= Hnnz+1 !convencao 


	  !matriz de precondicionamento DIAGONAL  !ccs, monta por coluna
	  Hnnz = 0 ! numero de elementos nao nulos da matriz de precondicionamento 
	  do jL=1,Nz0  !loop no grid de incognitas do pi
	      do iL=1,Nx0 
		  noL = (jL-1)*Nx0 + il !numero da incognita (coluna da Hessiana)

		  Hnnz=Hnnz+1  !H(L,L)
		  Hcolptr(noL)= Hnnz	   ! nnz do primeiro nao nulo da coluna ! CCS da H 
		  Hrowind(Hnnz)= noL   ! numero da linha! CCS da H 
		  Hcolind(Hnnz)= noL     ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
		  Hvalue(Hnnz)=0.0d0   ! CCS da H 
	      enddo
	  enddo
	  Hcolptr(Nx0*Nz0+1)= Hnnz+1 !convencao 

! ! 	  matriz de precondicionamento 3DIAGONAIS'  !ccs, monta por coluna
! 	  Hnnz = 0 ! numero de elementos nao nulos da matriz de precondicionamento
! 	  do jL=1,Nz0  !loop no grid de incognitas do pi
! 	      do iL=1,Nx0 
! 		  noL = (jL-1)*Nx0 + il !numero da incognita (coluna da Hessiana)
! 			
! 		  if(iL.ne.1 .and. iL.ne. Nx0 .and. jL.ne.1 .and. jL.ne.Nz0)then
! 
! 		      Hnnz=Hnnz+1 !H(L-1,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL ! numero da  linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 		 
! 		  elseif(iL.eq.1 .and. jl.eq.1)then
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(iL.eq.1 .and. jl.eq.Nz0)then
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(iL.eq.Nx0 .and. jl.eq.Nz0)then
! 
! 		      Hnnz=Hnnz+1 !H(L-1,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(iL.eq.Nx0 .and. jl.eq.1)then
! 
! 		      Hnnz=Hnnz+1  !H(L-1,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(iL.eq.1)then
! 
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
!  		  elseif(iL.eq.nX0)then
! 
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(jL.eq.1)then
! 
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(jL.eq.Nz0)then
! 
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 	
! 		  endif
! 
! 	      enddo
! 	  enddo ! fim do loop das incognitas
! 	  Hcolptr(Nx0*Nz0+1)= Hnnz+1 !convencao 

! ! 	  matriz de precondicionamento 5DIAGONAIS'  !ccs, monta por coluna
! 	  Hnnz = 0 ! numero de elementos nao nulos da matriz de precondicionamento
! 	  do jL=1,Nz0  !loop no grid de incognitas do pi
! 	      do iL=1,Nx0 
! 		  noL = (jL-1)*Nx0 + il !numero da incognita (coluna da Hessiana)
! 
! 		  if(iL.ne.1 .and. iL.ne. Nx0 .and. jL.ne.1 .and. jL.ne.Nz0)then
! 
! 		      Hnnz=Hnnz+1  !H(L-nx0,L)
! 		      Hcolptr(noL)= Hnnz	   ! nnz do primeiro nao nulo da coluna ! CCS da H 
! 		      Hrowind(Hnnz)= noL-Nx0   ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL     ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0   ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L-1,L)
! 		      Hrowind(Hnnz)= noL-1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL ! numero da  linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L+nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 	
! 
! 		  elseif(iL.eq.1 .and. jl.eq.1)then
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+Nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(iL.eq.1 .and. jl.eq.Nz0)then
! 		      Hnnz=Hnnz+1 !H(L-Nx0,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(iL.eq.Nx0 .and. jl.eq.Nz0)then
! 		      Hnnz=Hnnz+1 !H(L-Nx0,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1 !H(L-1,L)
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 
! 
! 		  elseif(iL.eq.Nx0 .and. jl.eq.1)then
! 		      Hnnz=Hnnz+1  !H(L-1,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1  !H(L+nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(iL.eq.1)then
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
!  
! 		  elseif(iL.eq.nX0)then
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(jL.eq.1)then
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(jL.eq.Nz0)then
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 	
! 		  endif
! 
! 	      enddo
! 	  enddo ! fim do loop das incognitas
! 	  Hcolptr(Nx0*Nz0+1)= Hnnz+1 !convencao 


! ! 	  matriz de precondicionamento 9DIAGONAIS'  !ccs, monta por coluna
! 	  Hnnz = 0 ! numero de elementos nao nulos da matriz de precondicionamento
! 	  do jL=1,Nz0  !loop no grid de incognitas do pi
! 	      do iL=1,Nx0 
! 		  noL = (jL-1)*Nx0 + il !numero da incognita (coluna da Hessiana)
! 
! 		  if(iL.ne.1 .and. iL.ne. Nx0 .and. jL.ne.1 .and. jL.ne.Nz0)then
! 
! 		      Hnnz=Hnnz+1  !H(L-nx0,L)
! 		      Hcolptr(noL)= Hnnz	   ! nnz do primeiro nao nulo da coluna ! CCS da H 
! 		      Hrowind(Hnnz)= noL-Nx0-1   ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL     ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0   ! CCS da H 
! 
! 		      Hnnz=Hnnz+1  !H(L-nx0,L)
! 		      Hrowind(Hnnz)= noL-Nx0   ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL     ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0   ! CCS da H 
! 
! 		      Hnnz=Hnnz+1  !H(L-nx0,L)
! 		      Hrowind(Hnnz)= noL-Nx0+1   ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL     ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0   ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L-1,L)
! 		      Hrowind(Hnnz)= noL-1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL ! numero da  linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 
! 
! 		      Hnnz=Hnnz+1 !H(L+nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0-1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 	
! 
! 		      Hnnz=Hnnz+1 !H(L+nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 	
! 
! 		      Hnnz=Hnnz+1 !H(L+nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0+1 ! numero da linha! CCS da H 
! 		      Hcolind(Hnnz)= noL ! numero da coluna ! esse eh para facilitar fazer o somatorio das derivadas
! 		      Hvalue(Hnnz)=0.0d0 ! CCS da H 	
! 
! 		  elseif(iL.eq.1 .and. jl.eq.1)then
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+Nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1 !H(L+Nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0+1
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(iL.eq.1 .and. jl.eq.Nz0)then
! 		      Hnnz=Hnnz+1 !H(L-Nx0,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1 !H(L-Nx0,L)
! 		      Hrowind(Hnnz)= noL-Nx0+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L+1,L)
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		  elseif(iL.eq.Nx0 .and. jl.eq.Nz0)then
! 
! 		      Hnnz=Hnnz+1 !H(L-Nx0,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0-1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1 !H(L-Nx0,L)
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1 !H(L-1,L)
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 
! 
! 		  elseif(iL.eq.Nx0 .and. jl.eq.1)then
! 		      Hnnz=Hnnz+1  !H(L-1,L)
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1 !H(L,L)
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1  !H(L+nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0-1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1  !H(L+nx0,L)
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(iL.eq.1)then
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-Nx0+1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
!  
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0+1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(iL.eq.nX0)then
! 
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0-1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0-1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(jL.eq.1)then
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0-1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+Nx0+1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		  elseif(jL.eq.Nz0)then
! 
! 		      Hnnz=Hnnz+1
! 		      Hcolptr(noL)= Hnnz
! 		      Hrowind(Hnnz)= noL-Nx0-1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-Nx0	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-Nx0+1	
! 		      Hcolind(Hnnz)= noL	  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL-1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0	
! 
! 		      Hnnz=Hnnz+1
! 		      Hrowind(Hnnz)= noL+1
! 		      Hcolind(Hnnz)= noL		  
! 		      Hvalue(Hnnz)=0.0d0
! 	
! 		  endif
! 
! 	      enddo
! 	  enddo ! fim do loop das incognitas
! 	  Hcolptr(Nx0*Nz0+1)= Hnnz+1 !convencao 
! 
! 
! 	  !matriz de precondicioanmento completa (so para problemas pequenos, testes) 
! 	  do iL=1,Nx0*Nz0
! 	      do jL=1,Nx0*Nz0
! 		  Hprec(il,jl) = 0.0d0
! 	      enddo
! 	  enddo !fim da inicializacao da metriz de precondicionamento cheia


	  write(*,*)'matriz H zerada'

      endif ! fim do iflagH=1: inicializacao da matriz de precondicionamento




      do igsr=1,gruposSR !loop nos grupos de fonte/receptor

	  !Monta o vetor independente para cada fonte (um elemento nao nulo em cada vetor)
	  do ks=1,nsource(igsr)
	      do i=1,Neq !zera os vetores referentes as fontes
		  f(i,ks) = (0.0d0,0.0d0)
	      enddo
	  enddo
	  do ks=1,Nsource(igsr)
	    ndonoh = (jsource(ks,igsr) -1)*Nx + isource(ks,igsr) !numero do noh da fonte ks
	    f(ndonoh,ks) = Sw(kw)/(rho(isource(ks,igsr),jsource(ks,igsr))*delta**2)  ! a divisao por rho é para que esta equacao fique de acordo com a equacao da modelagem no tempo (que é aquela Uxx+Uzz-1/c**2 Utt=S)
	    !write(*,*) 'fonte',f((jfonte -1)*Nx + ifonte)
	  enddo

! 	      ! Solver: eliminacao de  gauss do lapack
! 	      NRHS = Nsource
! 	      allocate(IPIV(Neq))
!       ! 	write(*,*)'#Resolvendo sistema........'
! 	      call CGBSV(Neq,KL,KU,NRHS,LwB,LDAB,IPIV,f,Neq,INFO)! matriz em banda
!       ! 	call CGESV(Neq,NRHS,Lw,Neq,IPIV,f,Neq,INFO)
! 	      write(*,*)'fim da sol do sistema........'

	  ! Solver: eliminacao de  gauss SuperLU
	  nrhs = nsource(igsr)
	  ldf = Neq
              
	  if(igsr.eq.1)then !	* First, factorize the matrix.(iopt=1) The factors are stored in *factors* handle.
	      iopt = 1 ! fatorar a primeira vez
	      call c_fortran_zgssv( iopt, neq, nnz, nrhs, values, rowind, colptr, &
			  & f, ldf, factors, info ) !CCS
	      if (info .eq. 0) then
! 		  write (*,*) 'Factorization succeeded'
	      else
		  write(*,*) 'INFO from factorization = ', info
		  stop
	      endif
	  endif


	  iopt = 2! aproveitar a matriz ja fatorada para resolver o sistema
	  call c_fortran_zgssv( iopt, neq, nnz, nrhs, values, rowind, colptr, &
			      & f, ldf, factors, info ) !CCS
	  if (info .eq. 0) then
! 	      write (*,*) 'Solve succeeded'
	  else
	      write(*,*) 'INFO from Solve = ', info
	      stop
	  endif


	  ! salva a resposta em todo o dominio, para todas as fontes, para a freq kw 
	  !(só serve para fazer snapshot)
	  ideltaSnap = int(Nsource(igsr)/(Nsourcesnap-1))
	  if (iflagF.eq.0 .and. iflagG.eq.0 .and. iflagH.eq.0)then! soh simulacao (snapshot)
	      do kss=1,NsourceSnap
		  ks=(kss-1)*ideltaSnap + 1
		  if(ks.gt.Nsource(igsr))ks=Nsource(igsr)
		  do i=1,Nx
		      do j=1,Nz
			  linha = (j-1)*Nx + i ! numero do noh i,j
			  P(i,j,kss,igsr) =  f(linha,ks)
		      enddo
		  enddo
	      enddo
	  endif


	  ! salva a resposta nos receptores, para todas as fontes, para a freq kw 
	  !(serve para fazer a inversao e p gerar sismograma)
	  do ks=1,Nsource(igsr)
	      do kr=1,Nrec(igsr)
		  linha = (jReceiv(kr,igsr)-1)*Nx + iReceiv(kr,igsr)  ! numero do noh do receptor i
		  Preceiv(kr,ks,kw,igsr) =  f(linha,ks)
	      enddo	
	  enddo

	
	  if(iflagF.eq.1)then

	      ! residuo complexo
	      do ks=1,Nsource(igsr)
		  do kr=1,Nrec(igsr)
		      resid((ks-1)*Nrec(igsr)+kr) = PreceivBARRA(kr,ks,kw,igsr) - Preceiv(kr,ks,kw,igsr)
		  enddo
	      enddo

	      do ks=1,Nsource(igsr)
		  do kr=1,Nrec(igsr)
		      Misfit = Misfit + resid((ks-1)*Nrec(igsr)+kr)*conjg(resid((ks-1)*Nrec(igsr)+kr))
		  enddo
	      enddo

	  endif

	  if(iflagG.eq.1)then

	      ! write(*,*)'calculo do grad da Misfit'
	      do ks=1,Nsource(igsr)
		  do i=1,Nx*Nz
		      residexpandido(i,ks) = 0.0 ! zera a matrix de vetores "residuo expandido" referente a acada fonte
		  enddo
	      enddo

	      do ks=1,Nsource(igsr)
		do kr=1,Nrec(igsr)
		      ndonoh = (jReceiv(kr,igsr)-1)*Nx + iReceiv(kr,igsr)  ! numero do noh do receptor kr
		      residexpandido(ndonoh,ks) =  CONJG(resid((ks-1)*Nrec(igsr)+kr))  ! residuo conjugado, (fonte adjunta)
		enddo
	      enddo

	      !NRHS = Nsource
	      !CALL CGBTRS('No transpose', Neq, KL, KU, NRHS, LwB, LDAB, IPIV, residexpandido, Neq, INFO )


	      !solve the system using the existing factors. SuperLU
	      nrhs = nsource(igsr)
	      ldf = Neq
	      iopt = 2
	      call c_fortran_zgssv( iopt, neq, nnz, nrhs, values, rowind, colptr, &
				    & residexpandido, ldf, factors, info ) !CCS
	      if (info .eq. 0) then
! 		  write (*,*) 'Solve succeeded _ grad'
	      else
		  write(*,*) 'INFO from triangular solve (p/ cal grad) = ', info
		  stop
	      endif

  
	      do noL=1,Nx0*nz0 ! loop nas variaveis do PI

		  cont = 1
		  do while(depX(cont,noL).ne.0)
		      ndonoh = depX(cont,noL)
		      i = mod(ndonoh,Nx)
		      if(i.eq.0)i=Nx
		      j= (ndonoh-1)/Nx + 1

		      ! agora com mais fontes:
		      aux = 0.0d0
		      do ks=1,Nsource(igsr)
			  aux = aux + REAL( f(ndonoh,ks) *  residexpandido(ndonoh,ks))
		      enddo
		      gradmisfit(noL) = gradmisfit(noL) + &
				  & (2.0d0*(Cmax-c(i,j))*(c(i,j)-Cmin)/(Cmax-Cmin))* &
				  &  (-2.0d0*omega**2/(rho(i,j)*c(i,j)**3)) *  aux      
		      cont=cont+1
		  enddo

	      enddo ! fim do loop nas variaveis do PI 
! 	      write(*,*)' grad form adjunta: ', gradmisfit  
! 	      fim do calculo do grad Misfit' ---- no fim do loop igsr ainda tem que dividir o vetor gradiete pelo denominador constante
	  endif !fim do if iflagG.eq.1


	  if(iflagH.eq.1)then

	      do kr=1,Nrec(igsr) ! zera a matrix du/dc referente a acada fonte
		  do i=1,nx*nz
		      deltaR(i,kr) = 0.0d0 ! zera todos os elemntos dos vetores coluna
		  enddo
		  ndonoh = (jReceiv(kr,igsr)-1)*Nx + iReceiv(kr,igsr)  ! numero do noh do receptor kr
		  deltaR(ndonoh,kr)= 1.0d0 !coloca 1.0 no noh do receptor (cada coluna tem um unico elemento diferente de zero)
	      enddo
	      write(*,*)'matriz deltaR criada'

	      write(*,*)'backpropagar deltaR:'
	      !solve the system using the existing factors. SuperLU
	      nrhs = Nrec(igsr)
	      ldf = Neq
	      iopt = 2
	      call c_fortran_zgssv( iopt, neq, nnz, nrhs, values, rowind, colptr, &
				    & deltaR, ldf, factors, info ) !CCS
	      if (info .eq. 0) then
! 		  write (*,*) 'Solve succeeded _ deltaR'
	      else
		  write(*,*) 'INFO from triangular solve (p/ calc de J) = ', info
		  stop
	      endif


	      do ks=1,Nsource(igsr)

		  do noL=1,Nx0*nz0 ! loop no numero de incognitas do problema inverso

		      do kr=1,Nrec(igsr)
			JdeX(kr,noL) = 0.0d0 ! esta matriz jacobiana eh gerada para cada grupo de pares fonte receptor
		      enddo

		      cont = 1
		      do while(depX(cont,noL).ne.0)
			  ndonoh = depX(cont,noL)
			  i = mod(ndonoh,Nx)
			  if(i.eq.0)i=Nx
			  j= (ndonoh-1)/Nx + 1

			  do kr=1,Nrec(igsr)
			      JdeX(kr,noL) = JdeX(kr,noL) + f(ndonoh,ks)* &
							&  (2.0d0*(Cmax-c(i,j))*(c(i,j)-Cmin)/(Cmax-Cmin))*&
							&  (-2.0d0*omega**2/(rho(i,j)*c(i,j)**3)) * &
							&  deltaR(ndonoh,kr) 
			  enddo
			  cont=cont+1
		      enddo !fim do do while

		  enddo ! fim do loop nas incognitas L

		    
		  do knnz=1,Hnnz ! matriz de precondicionamento H (so elementos da diagonal de J^t * J)
		      do kr=1,Nrec(igsr)
			  Hvalue(knnz)= Hvalue(knnz)+ DREAL(CONJG(JdeX(kr,Hcolind(knnz)))*JdeX(kr,Hrowind(knnz)))
		      enddo
! 			  write(1,*)Hrowind(knnz),Hcolind(knnz),Hvalue(knnz)
		  enddo !fim da matriz de precondicionamento'
		
! 		  do jL=1,Nx0*Nz0 ! montagem na matriz de precondicionamento completa J^T * J
! 		      do iL=jl,Nx0*Nz0
! 			  do kr=1,Nrec(igsr)
! 			      Hprec(il,jl) = Hprec(il,jl)+ REAL(CONJG(JtdeX(il,kr))*JtdeX(jl,kr))
! 			  enddo
! 		      enddo
! 		  enddo

	      enddo ! fim do loop das fontes

	  endif ! fim do if iflagH =1 hessiana

      enddo  !loop nos grupos de fonte/receptor




      if(iflagF.eq.1)then
! 	  Misfit = Misfit! funcao objetivo tradicional  sem multiplicar a norma^2 por 1/2
! 	  errorelat = dsqrt(Misfit/norm2Pbarra)
! 	  Misfit = Misfit/2.0d0 ! funcao objetivo tradicional
! 	  errorelat = dsqrt(Misfit*2.0d0/norm2Pbarra)
! 	  Misfit = Misfit/(2.0d0*norm2Pbarra) ! funcao objetivo normalizada sem ponderacao de hu2011 , eq. 8, sem a parcela da regularizacao
! 	  errorelat = dsqrt(Misfit*2.0d0)
	  Misfit = Misfit/(norm2Pbarra) ! funcao objetivo normalizada sem ponderacao de hu2011 , eq. 8, sem a parcela da regularizacao
	  errorelat = dsqrt(Misfit)
! 	  Misfit = Misfit/(2.0d0*medidas)
!  	  errorelat = dsqrt(Misfit*2.0d0*medidas/norm2Pbarra)
      endif
      if(iflagG.eq.1)then
	  do noL=1,Nx0*Nz0
! 	      gradmisfit(noL) = 2.0d0*gradmisfit(noL)! funcao objetivo tradicional  sem multiplicar a norma^2 por 1/2
! 	      gradmisfit(noL) = gradmisfit(noL)! grad da funcao objetivo tradicional
! 	      gradmisfit(noL) = gradmisfit(noL)/norm2Pbarra! grad da funcao objetivo normalizada sem ponderacao de hu2011
	      gradmisfit(noL) = 2.0d0*gradmisfit(noL)/norm2Pbarra! grad da funcao objetivo normalizada sem ponderacao de hu2011
! 	      gradmisfit(noL) = gradmisfit(noL)/dfloat(medidas)
	  enddo
      endif
! ! ! !       if(iflagH.eq.1)then
! ! ! ! 	  do noL=1,Nx0*Nz0
! ! ! ! 	      do noK=1,Nx0*Nz0
! ! ! ! 		  Hprec(noK,noL) = 0.0d0
! ! ! ! 	      enddo
! ! ! ! 	  enddo
! ! ! ! 	  do knnz=1,Hnnz ! matriz de precondicionamento H (so elementos da diagonal de J^t * J)
! ! ! ! 	      Hprec(Hrowind(knnz),Hcolind(knnz)) = (Hvalue(knnz))
! ! ! ! 	  enddo
! ! ! ! 	  call system('rm hessiana.bin')
! ! ! ! 	  open(1,file='hessiana.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx0*Nz0*Nx0*Nz0)
! ! ! ! 	  write(1,rec=1)(((Hprec(i,j)),i=1,nx0*nz0),j=1,nx0*nz0)
! ! ! ! 	  close(1)
! ! ! ! ! 	  stop
! ! ! ! 
! ! ! ! ! 	  !uma vez que H ja foi atualizada, é fatorada aqui:
! ! ! ! ! 	  nrhsH = 1
! ! ! ! ! 	  ldyH = Nx0*nz0
! ! ! ! ! 	  ioptH = 1 ! fatorar 
! ! ! ! ! 	  call c_fortran_dgssv( ioptH, ldyH, Hnnz, nrhsH, Hvalue, Hrowind, Hcolptr, &
! ! ! ! ! 		      & f, ldyH, factorsH, infoH ) !CCS !aqui o f nao é usado, esta rotina nao conhece o vetor independente do sistema, coloquuei qq um
! ! ! ! ! 	  if (infoH .eq. 0) then
! ! ! ! ! 	      write (*,*) 'H Factorization succeeded'
! ! ! ! ! 	  else
! ! ! ! ! 	      write(*,*) 'H INFO from factorization = ', infoH
! ! ! ! ! 	      stop
! ! ! ! ! 	  endif
! ! ! !       endif ! fim do if iflagH.eq.1




      ! * Last, free the storage allocated inside SuperLU
      iopt = 3
      call c_fortran_zgssv( iopt, neq, nnz, nrhs, values, rowind, colptr, &
				  & f, ldf, factors, info ) !CCS
      if (info .eq. 0) then
! 	  write (*,*) 'free the storage succeeded'
      else
	  write(*,*) 'INFO from free the storage = ', info
	  stop
      endif


!       ! lu do lapack
!       deallocate(IPIV)
!       deallocate(LwB)
!       deallocate(f)

      end subroutine PDMDFDF1